Projekt zakładał analizę danych pacjentów chorych na COVID-19. Przetworzenie danych ze zbioru umożliwiło sporządzenie raportu zawierającego wykresy takie jak rozkład wyzdrowień i śmierci. Podjęto także próbę stworzenia klasyfikatora, który będzie potrafił wskazać czy dany pacjent przeżyje.
library("readxl")
library("dplyr")
library(skimr)
library(plotly)
library(ggplot2)
library(DT)
library(hablar)
library(reshape2)
library(tidyr)
library(caret)
set.seed(23)Data próbki została przycięte do formatu “yyyy-MM-dd”.Następnie dane zostały zgrupowane pod dacie przyjęcia oraz dacie opuszenia szpitala, z powodu wyzdrowienia lub zgonu. Co umożliwiło przypisanie brakujących PATIENT_ID.
df <- mutate(data, RE_DATE = as.Date(substr(RE_DATE,1,10), "%Y-%m-%d"))
df <- df %>% group_by(Admission.time, Discharge.time) %>% mutate(PATIENT_ID = max(PATIENT_ID, na.rm = TRUE)) %>% ungroup()
df <- df %>% group_by(PATIENT_ID, RE_DATE) %>%
mutate(Admission.time = max(Admission.time, na.rm = TRUE)) %>%
mutate(Discharge.time = max(Discharge.time, na.rm = TRUE))%>% ungroup()Wszystkie próbki z danego dnia dla danego pacjenta zostały zgrupowane.
df <- df %>% group_by(PATIENT_ID, RE_DATE, Admission.time, Discharge.time, age, outcome, gender) %>% summarise_all(sum_) %>% ungroup()Dane liczbowe odpowiadające płci zostały zamienione na dane tekstowe.
df <- df %>% mutate(gender = factor(gender, labels = c('male', 'female')))
df <- df %>% mutate(outcome = factor(outcome, labels = c('survived', 'died')))Kolumny i wiersze zawierające niekompletne dane zostały usunięte.
df$X2019.nCoV.nucleic.acid.detection = NULL
df <- subset(df, !is.na(RE_DATE))Tabela przedstawia zbiór po przetworzeniu.
datatable(df,options = list(scrollX = TRUE), filter = 'top')my_skim <- skim_with(character = sfl(), append = FALSE, base = sfl(n_missing))
datatable(my_skim(df), options = list(scrollX = TRUE))Wykresy przedstawiają wpływ wartości dehydrogenazy mleczanowej (LDH), białka C-reaktywnego (hs-CRP), oraz liczby limfocytów.
last_data <- df
last_data <- last_data %>% group_by(PATIENT_ID) %>%
summarise(outcome = last(outcome), LDH = last(Lactate.dehydrogenase), Limfocyty = last(Lactate.dehydrogenase), 'hs-CRP' = last(High.sensitivity.C.reactive.protein)) %>% ungroup()
last_data[ ,c('PATIENT_ID')] <- list(NULL)
datatable(last_data)last_plot_data = melt(last_data, id=c("outcome"))
graph <- ggplot(last_plot_data, aes(x=value, y=outcome)) + geom_point() + facet_wrap(~variable)
ggplotly(graph)Wykres przedstawia dane dla każdego dnia. Zsumowana została liczba pacjentów opuszczających szpital z podziałem na zgony oraz wyzdrowienia.
by_date <- mutate(df, Discharge.time = substr(Discharge.time,1,10))
by_date <- by_date %>% group_by(Discharge.time, PATIENT_ID) %>% mutate(outcome = first(outcome == "died")) %>%
summarise(outcome = first(outcome)) %>% group_by(Discharge.time) %>%
summarise(Pacjenci = n_distinct(PATIENT_ID), Zgony = sum(outcome == TRUE), Wyzdrowienia = sum(outcome == FALSE))
plot_data = melt(by_date, id=c("Discharge.time"))
p <- ggplot(plot_data) +
geom_point(aes(x = Discharge.time, y = value, color = variable,
text = paste(' Dzień: ', Discharge.time,'<br>','Liczba: ',value,'<br>','Typ: ',variable))) +
scale_colour_manual(values=c("black","red","green")) +
labs(color = "Legenda", x = "Dzień", y = "Liczba") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
g <- ggplotly(p, tooltip = "text") %>% layout(legend = list(x = 0.05, y = 0.95))
gPrzed przystąpienie do uczenia klasyfikatora dane zostały odpowiednio przygotowane. Ze zbioru zostali usunięci pacjenci ze zbyt mała liczbą przeprowadzonych badan oraz badańia wykonane na zbyt małej liczbie pacjentów.
raw_train_data <- df
raw_train_data[ ,c('Admission.time', 'Discharge.time')] <- list(NULL)
raw_train_data <- raw_train_data[which(rowMeans(!is.na(raw_train_data)) > 0.7), which(colMeans(!is.na(raw_train_data)) > 0.32)]Po wstępnym przetworzeniu zbioru, usunięto wszystkich pacjentów, którzy nie mieli kompletu pozostałych badań.
raw_train_data <- raw_train_data %>% drop_na()
patient_data <- raw_train_data %>% group_by(PATIENT_ID) %>% summarise_at(vars(outcome, gender), last) %>% ungroup()
raw_train_data <- raw_train_data %>% group_by(PATIENT_ID) %>% summarise_if(is.numeric, max) %>% ungroup()Tabela przedstawia ostateczny zbiór, który zostanie wykorzystany do uczenia klasyfikatora. Zawiera on 217 rekordóW.
raw_train_data <- merge(x = patient_data, y = raw_train_data, by ="PATIENT_ID", all = TRUE)
datatable(raw_train_data,options = list(scrollX = TRUE))inTraining <- createDataPartition(y = raw_train_data$outcome, p = .75, list = FALSE)
training <- raw_train_data[inTraining,]
testing_val <- raw_train_data[-inTraining,]
testing_val_index <- createDataPartition(y = testing_val$outcome, p = .6, list = FALSE)
testing <- testing_val[-testing_val_index,]
val <- testing_val[testing_val_index,]
datatable(training, options = list(scrollX = TRUE))datatable(testing, options = list(scrollX = TRUE))datatable(val, options = list(scrollX = TRUE))Dane zostały wykorzystane do uczenia klasyfikator random forest.
ctrl <- trainControl(method = "repeatedcv", number = 2,repeats = 5)
fit <- train(outcome ~ .,
data = training,
method = "rf",
trControl = ctrl,
ntree = 10)
fit## Random Forest
##
## 163 samples
## 58 predictor
## 2 classes: 'survived', 'died'
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 81, 82, 81, 82, 82, 81, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.9164408 0.8305415
## 30 0.9839958 0.9676142
## 58 0.9950918 0.9901007
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 58.rfGrid <- expand.grid(mtry = 10:30)
gridCtrl <- trainControl(method = "repeatedcv", summaryFunction = twoClassSummary, classProbs = TRUE, number = 2, repeats = 5)
fitTune <- train(outcome ~ .,
data = training,
method = "rf",
metric = "ROC",
preProc = c("center", "scale"),
trControl = gridCtrl,
tuneGrid = rfGrid,
ntree = 30)
fitTune## Random Forest
##
## 163 samples
## 58 predictor
## 2 classes: 'survived', 'died'
##
## Pre-processing: centered (58), scaled (58)
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 81, 82, 82, 81, 81, 82, ...
## Resampling results across tuning parameters:
##
## mtry ROC Sens Spec
## 10 0.9933192 0.9666667 0.9669670
## 11 0.9986728 0.9777778 0.9835586
## 12 0.9976852 0.9688889 0.9807057
## 13 0.9978095 0.9755556 0.9863363
## 14 0.9971330 0.9688889 0.9725976
## 15 0.9971296 0.9755556 0.9694444
## 16 0.9979630 0.9844444 0.9780030
## 17 0.9983025 0.9755556 0.9888889
## 18 0.9981790 0.9666667 1.0000000
## 19 0.9982716 0.9777778 0.9917417
## 20 0.9991358 0.9822222 0.9945195
## 21 0.9984259 0.9733333 0.9833333
## 22 0.9989506 0.9888889 0.9889640
## 23 0.9983959 0.9844444 0.9972973
## 24 0.9988889 0.9822222 0.9889640
## 25 0.9992009 0.9866667 1.0000000
## 26 0.9990741 0.9844444 0.9972222
## 27 0.9999074 0.9866667 1.0000000
## 28 0.9991975 0.9844444 1.0000000
## 29 0.9996914 0.9888889 0.9972973
## 30 0.9992593 0.9888889 1.0000000
##
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 27.Uzyskano bardzo wysoki wskaźnik “accuracy”. Nauczony klasyfikator poprawnie sklasyfikował wszystkie próbki. Niestety liczba danych wykorzystana w procesie jest za mała. Nie można mieć pewności, że klasyfikator będzie działać z tak dobrą skutecznością na znacznie większym zbiorze.
predict_result <- predict(fitTune, newdata = val)
confusionMatrix(data = predict_result, val$outcome)## Confusion Matrix and Statistics
##
## Reference
## Prediction survived died
## survived 18 0
## died 0 15
##
## Accuracy : 1
## 95% CI : (0.8942, 1)
## No Information Rate : 0.5455
## P-Value [Acc > NIR] : 2.056e-09
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0000
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 1.0000
## Prevalence : 0.5455
## Detection Rate : 0.5455
## Detection Prevalence : 0.5455
## Balanced Accuracy : 1.0000
##
## 'Positive' Class : survived
##